#include "lineardiscriminantanalysis.h"

LinearDiscriminantAnalysis::LinearDiscriminantAnalysis(int feat)
{
    src = 0;
    projMat = 0;
    class_info = new vector<int>;
    numClass = 2;
    f_rows = feat;
    d_rows = 0;
}

LinearDiscriminantAnalysis::~LinearDiscriminantAnalysis()
{
    if(src) cvReleaseMat(&src);
    if(projMat) cvReleaseMat(&projMat);
    projMat = 0;
    src = 0;

    if(class_info) class_info->clear();
    delete class_info;
}

/*
 *  Read file and store data into cvMat* src
 *
 */
bool LinearDiscriminantAnalysis::loaddata(const char *file)
{
    //open file for reading
    fstream inputFile;
    inputFile.open(file, ios::in);

    if ( inputFile.is_open() )
    {
        // Temporary holding data from input file
        vector<double *> temp;

        string line = "";

        //read and store data
        while ( !inputFile.eof() )
        {
            getline(inputFile, line);

            //process line
            if (line.length() > 2 )
            {
                processLine(temp, line);
                // number of records/lines
                d_rows++;
            }
        }

        // Store data into cvMat
        storeData(temp, d_rows);

        // Reduce dimension of feature vector
        projectionMat();

        // Output projectionMat
        writeData(projMat);

        // clear memory for input vector
        for(int i = 0; i < d_rows; i++)
            delete [] temp.at(i);
        temp.clear();

        inputFile.close();
        return true;
    }
    return false;
}

// Process input line by line from file and save each line into vector
void LinearDiscriminantAnalysis::processLine(vector<double *> &temp, string &line)
{
    //create new pattern
    double* pattern = new double[f_rows+1];

    //store inputs
    char* cstr = new char[line.size()+1];
    char* t;
    strcpy(cstr, line.c_str());

    //tokenise
    int i = 0;
    t=strtok (cstr,",");

    while ( t!=NULL && i < (f_rows+1) )
    {
        pattern[i++] = atof(t);

        //move token onwards
        t = strtok(NULL,",");
    }

    //add to records
    temp.push_back( pattern );
    delete [] cstr;
}
// After getting input from input file, then store them into matrix 'src'
void LinearDiscriminantAnalysis::storeData(vector<double *> &temp, int rows)
{
    if(rows <= 0) return;

    int c = 0;
    // srcmat to hold input values in rotated dimension
    src = cvCreateMat(f_rows, rows, CV_64FC1);

    for(int i = 0; i < rows; i++)
    {
        double *line = temp.at(i);
        int j = 0;
        for(; j < f_rows; j++)
        {
            cvmSet(src, j, i, line[j]);
        }
        c = ((int)line[j] > 0) ? 1 : 0;
        class_info->push_back(c);
    }
}

void LinearDiscriminantAnalysis::projectionMat()
{
    CvMat *scatter_within = cvCreateMat(f_rows, f_rows, CV_64FC1);
    CvMat* scatter_between = cvCreateMat(f_rows, f_rows, CV_64FC1);
    CvMat* classMean = cvCreateMat(f_rows, numClass, CV_64FC1);

    classMeanMat(classMean);
    withinMat(classMean, scatter_within);
    betweenMat(classMean, scatter_between);

    // Projection Mat
    ratioMat(scatter_between, scatter_within);

    cvReleaseMat(&scatter_within);
    cvReleaseMat(&scatter_between);
}

// Calculate mean within each class
void LinearDiscriminantAnalysis::withinMat(CvMat *classMean, CvMat *within)
{
    int index;
    double temp = 0.0;

    CvMat* correctedMean = cvCreateMat(f_rows, d_rows, CV_64FC1);

    for(int i = 0; i < d_rows; i++)
    {
        index = class_info->at(i);
        for(int j=0; j<f_rows; j++)
        {
            temp = cvmGet(src, j, i) - cvmGet(classMean, j, index);
            cvmSet(correctedMean, j, i, temp);
        }
    }

    cvGEMM(correctedMean, correctedMean, 1.0, 0, 0.0, within, CV_GEMM_B_T);

    cvReleaseMat(&correctedMean);
}

void LinearDiscriminantAnalysis::betweenMat(CvMat *classMean, CvMat *between)
{
    covarianceMat(classMean, between);
}

void LinearDiscriminantAnalysis::ratioMat(CvMat *s_b, CvMat *s_w)
{
    int i, j, index=0;
    double temp=0;

    CvMat* ratioMat = cvCreateMat(s_w->cols, s_b->cols, CV_64FC1);
    CvMat* inv_s_w = cvCreateMat(s_w->cols, s_w->rows, CV_64FC1);
    CvMat* W = cvCreateMat(s_w->cols, 1, CV_64FC1);
    CvMat* U = cvCreateMat(s_w->cols, s_w->cols, CV_64FC1);

    cvInvert(s_w, inv_s_w, CV_SVD);
    cvGEMM(inv_s_w, s_b, 1.0, 0, 0.0, ratioMat, 0);

    cvSVD(ratioMat, W, U, NULL, 0);

    for(i = 0; i <s_w->cols; i++)
    {
        if(cvmGet(W, i, 0) == 0.0)
        {
            index = i;
            break;
        }
    }

    projMat = cvCreateMat(s_w->cols, index+1, CV_64FC1);
    for(i=0; i<s_w->cols; i++)
    {
        for(j=0; j<index+1; j++)
        {
            temp = cvmGet(U, i, j);
            cvmSet(projMat, i, j, temp);
        }
    }

    cvReleaseMat(&ratioMat);
    cvReleaseMat(&inv_s_w);
    cvReleaseMat(&W);
    cvReleaseMat(&U);
}

void LinearDiscriminantAnalysis::classMeanMat(CvMat *mean)
{
    int count = 0;
    double tmp_val = 0.0;

    CvMat* sum = cvCreateMat(f_rows, 1, CV_64FC1);
    CvMat* col = cvCreateMat(f_rows, 1, CV_64FC1);

    for(int i=0;i<numClass;i++)
    {
        cvSetZero(sum);
        count=0;

        for(int j=0;j<d_rows;j++)
        {
          if(class_info->at(j) == i)
            {
              count++;
              cvGetCol(src,col,j);
              cvAdd(col,sum,sum,0);
            }
        }

        for(int j=0;j<f_rows;j++)
        {
           tmp_val=cvmGet(sum,j,0);
           cvmSet(mean,j,i,tmp_val/(double)count);
        }
    }
    cvReleaseMat(&sum);
    cvReleaseMat(&col);
}

void LinearDiscriminantAnalysis::covarianceMat(CvMat *classMean, CvMat *scatterMat)
{
    CvMat* mean_sub = cvCreateMat(f_rows, numClass, CV_64FC1);
    CvMat* mean = cvCreateMat(f_rows, 1, CV_64FC1);

    totalMeanMat(classMean, mean);

    int i, j;
    double temp;

    for(i = 0; i<numClass; i++)
    {
        for(j = 0; j<f_rows; j++)
        {
            temp = cvmGet(classMean, j, i) - cvmGet(mean, j, 0);
            cvmSet(mean_sub, j, i, temp);
        }
    }

    cvGEMM(mean_sub, mean_sub, 1.0, 0, 0.0, scatterMat, CV_GEMM_B_T);

    cvReleaseMat(&mean_sub);
    cvReleaseMat(&mean);
}

void LinearDiscriminantAnalysis::totalMeanMat(CvMat *classMean, CvMat *totalMean)
{
    int i, j;
    double sum = 0.0;

    for(i = 0; i < f_rows; i++)
    {
        sum = 0.0;
        for(j = 0; j < numClass; j++)
        {
            sum += cvmGet(classMean, i, j);
        }
        cvmSet(totalMean, i, 0, (sum/(double)numClass));
    }
}

bool LinearDiscriminantAnalysis::writeData(CvMat *data)
{
    if(!data) return false;

    int rows = data->rows;
    int cols = data->cols;

    fstream outputFile;
    if(!outputFile.is_open()) outputFile.open("./output/lda.csv", ios::out);

    if(outputFile.is_open())
    {
        double temp = 0;
        for(int i=0;i<rows;i++)
        {
            for(int j=0;j<cols;j++)
            {
               temp=cvmGet(data,i,j);
               outputFile << temp << ",";
            }
            outputFile << endl;
        }

        outputFile.close();
        return true;
    }
    return false;
}
